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Abstract 

We present a new Monte Carlo method for obtaining solutions of the Boltzmann 
equation for describing phonon transport in micro and nanoscale devices. The pro- 
posed method can resolve arbitrarily small signals (e.g. temperature differences) at 
small constant cost and thus represents a considerable improvement compared to 
traditional Monte Carlo methods whose cost increases quadratically with decreas- 
ing signal. This is achieved via a control-variate variance reduction formulation in 
which the stochastic particle description only solves for the deviation from a nearby 
equilibrium, while the latter is described analytically. We also show that simulating 
an energy-based Boltzmann equation results in an algorithm that lends itself nat- 
urally to exact energy conservation thereby considerably improving the simulation 
fidelity. Simulations using the proposed method are used to investigate the effect 
of porosity on the effective thermal conductivity of silicon. We also present simula- 
tions of a recently developed thermal conductivity spectroscopy process. The latter 
simulations demonstrate how the computational gains introduced by the proposed 
method enable the simulation of otherwise intractable multiscale phenomena. 



1 Introduction 

Over the past two decades, the dramatic advances associated with MEMS (Micro Elec- 
tro Mechanical Systems) and NEMS (Nano Electro Mechanical Systems) have attracted 
considerable attention on microscale and nanoscale heat transfer considerations [lj. Ap- 
plications range from thermal management of electronic devices [2] to the development 
of thermoelectric materials with higher figure of merit [3] . The thermoelectric figure of 
merit is proportional to the electrical conductivity and inversely proportional to ther- 
mal conductivity and can thus be improved by reducing the latter and/or increasing 
the former. One of the most promising approaches towards reducing the thermal con- 
ductivity of thermoelectric materials is the introduction of nanostructures that interact 



with the ballistic motion of phonons at small scales thus influencing heat transport |4J. 
Such approach requires a reliable description of phonon transport at the nanoscale and 
cannot rely on Fourier's law, which is valid for diffuse transport. On the other hand, 
first principles calculations (e.g. molecular dynamics approaches, classical or quantum 
mechanical) are too expensive for treating phonon transport at the device (e.g. microm- 
eter) scale. At these scales, a kinetic description based on the Boltzmann Transport 
Equation (BTE) offers a reasonable balance between fidelity and model complexity and 
is able to accurately describe the transition from diffusive to ballistic transport as char- 
acteristic system lengthscales approach and ultimately become smaller than the phonon 
mean free path. 

Solving the BTE is a challenging task, especially in complex geometries. The high 
dimensionality of the distribution function coupled with the ability of particle methods 
to naturally simulate advection processes without stability problems [5] make particle 
Monte Carlo methods particularly appealing. Following the development of the Direct 
Monte Carlo Method by Bird [6] for treating dilute gases, Monte Carlo methods for 
phonon transport were first introduced by Peterson [7] and subsequently improved by 
Mazumder and Majumdar [8]. Over the past decade, further important refinements have 
been introduced: Lacroix et al. introduced a method to treat frequency dependent mean 
free paths [9]; Jeng et al. introduced a method for efficiently treating transmission and 
reflection of phonons at material interfaces and used this method to model the thermal 
conductivity of nanoparticle composites [3J. Hao et al. developed |10j a formulation 
for periodic boundary conditions in order to study the thermal conductivity of periodic 
nanoporous materials while only simulating one unit cell (period). 

The work presented here introduces a number of improvements which enable ef- 
ficient and accurate simulation of the most challenging phonon transport problems, 
namely three-dimensional and transient. Accuracy is improved compared to previous 
approaches by introducing an energy-based formulation, which simulates energy packets 
rather than phonons; this formulation makes energy conservation particularly easy to 
implement rigorously, in contrast to previous approaches which were ad-hoc and in many 
cases ineffective. We also introduce a variance-reduced formulation for substantially re- 
ducing the statistical uncertainty associated with sampling solution (temperature and 
heat flux) fields. This formulation is based on the concept of control variates, first in- 
troduced in the context of Monte Carlo solutions of the Boltzmann equation for dilute 
gases p$j; it is based on the fact that signal strength is intimately linked to deviation 
from equilibrium, or in other words, that the large computational cost associated with 
small signals is due to the fact that in these problems the deviation from equilibrium 
is small. This observation can be exploited by utilizing the nearby equilibrium state 
as a "control" and using the Monte Carlo method to calculate the contribution of non- 
equilibrium therefrom. Because the deviation from equilibrium is small, only a small 
quantity is evaluated stochastically (the fields associated with the equilibrium compo- 
nent are known analytically) resulting in small statistical uncertainty; moreover, the 
latter decreases as the deviation from equilibrium decreases, thus enabling the simula- 
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tion of arbitrarily small deviations from equilibrium. 

In the technique presented here, we use particles to simulate the deviation from equi- 
librium, which is thus referred to as a deviational method; the origin of this methodol- 
ogy can be found in the Low Variance Deviational Simulation Monte Carlo (LVDSMC) 
method [HJ El US E] recently developed for dilute gases. The theoretical basis un- 
derlying this method as well as the modifications required for use in phonon transport 



simulations are described in section 2.3 The resulting algorithm is described in section 
[3] and validated in section HI 

The proposed algorithm is used to obtain solutions to two problems of practical inter- 
est. The first application studies the thermal conductivity of porous silicon containing 
voids with different degrees of alignment and is intended to showcase how ballistic ef- 
fects influence the "effective" thermal conductivity. The second application is related 
to the recently developed experimental method of "thermal conductivity spectroscopy" 
|15j based on the pump-probe technique known as transient thermoreflectance, which 
uses the response of a material to laser irradiation to infer information about physical 
properties of interest [16] (e.g. mean free paths of the dominant heat carriers). 

2 Theoretical basis 

2.1 Summary of traditional Monte Carlo simulation methods 

We consider the Boltzmann Transport Equation in the frequency-dependent relaxation- 
time approximation 

d f f - f loc 

where, / = f(t, r,k,p) is the phonon distribution in the phase space, oj = oj(\i,p) the 
phonon radial frequency, p the phonon polarization and T the temperature; similarly to 
the nomenclature adopted in p], / is defined in reference to the occupation number. For 
example, if the system is perfectly thermalized at temperature T eq , f is a Bose-Einstein 
distribution 

f e <? - (0) 

/t "-p(^)-i " 

where kb is Boltzmann's constant. Also, f loc is an equilibrium (Bose-Einstein) distri- 
bution corresponding to the local pseudo-temperature defined more precisely in section 
EX2J 

In this work we consider Longitudinal Acoustic (LA), Transverse Acoustic (TA), 
Longitudinal Optical (LO), Transverse Optical (TO) polarizations; acoustic phonons are 
known to be the most important contributors to lattice thermal conductivity \17\ 118] . 
The phonon radial frequency is given by the dispersion relation oj = oj(k,p). Phonons 
travel at the group velocity V g = V^w. 

In the following, we always consider the ideal case where the dispersion relation 
is isotropic. For convenience, the radial frequency oj and two polar angles 9 and cf) 
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are usually preferred as primary parameters compared to the wave vector. Equation 
([!]) is simulated using computational particles that represent phonon bundles, namely 
collections of phonons with similar characteristics (position vector x, the wave vector 
k, and the polarization/propagation- mode p), using the approximation 

±f(t, x, k,p) « N eff 5 3 (x - Xi )<5 3 (k - k*)^ (3) 

i 

where Xj, kj and pi respectively represent the position, the wave vector and the polar- 
ization of particle i and N e ff is the number of phonons in each phonon bundle. The 
factor 1/87T 3 is necessary for converting the quantity representing the occupation num- 
ber, /, into a quantity representing the phonon density in phase space. Written in 
polar coordinates, and using the frequency instead of the wave number, this expression 
becomes 



4ir 



f(t, x, W) 9, (t>,p) sin(0) « N eff <5 3 (x - ^)5(u - ^5(9 - 0<)«5(0 - 4>i)8 PjPi (4) 



where Ui, 6i, and fa respectively represent the radial frequency, the polar angle and the 
azimuthal angle of particle i. The density of states, D(oj,p), is made necessary by the 
use of uj as a primary parameter and is given by 

2.1.1 Initialization 

Systems are typically initialized in an equilibrium state at temperature T; the number 
of phonons in a given volume V is calculated using the Bose-Einstein statistics 

P^max 

N = V ^D(cj,p)f^(u)dw (6) 

J u;=0 p 

where: 

• u m ax is the maximum (cut-off) frequency 

• is the occupation number at equilibrium at temperature T 

The number of computational particles (each representing a phonon bundle) is given 
by N/N e ff. The value of N e ff is determined by balancing computational cost (includ- 
ing storage) with the need for a sufficiently large number of particles for statistically 
meaningful results. 
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2.1.2 Time integration 

Once the system is initialized, the simulation proceeds by applying a splitting algorithm 
with timestep At. Integration for one timestep comprises of three substeps: 

• The advection substep during which bundle i moves by V fl jAt. 



The sampling substep during which the temperature (T) and pseudo-temperature 
(T[ oc ) are locally measured. They are calculated by inverting the local energy (E) 
and pseudo-energy (E) [10] relations 

t ■ / "-° p exp (,§.)- 1 

and 

E = N V ^ = y r™*^ D{u,p)hw 1 
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respectively. 

• The scattering substep, during which each phonon i is scattered according to its 
scattering probability given by 

Pi = l-expf- At ) (9) 

Scattering proceeds by drawing new frequencies, polarizations and traveling di- 
rections. Because of the frequency dependent relaxation times, frequencies must 
be drawn from the distribution D(uj,p)f loc /T(uj,p,T). Since scattering events 
conserve energy, the latter must be conserved during this substep. However, be- 
cause the frequencies of the scattered phonons are drawn randomly, conservation 
of energy is enforced by adding or deleting particles until a target energy is ap- 
proximately reached [SI E]- In addition to being approximate, this method does 
not always ensure that energy is conserved, resulting in random walks in the en- 
ergy of the simulated system, which in some cases leads to deterministic error. In 
the next section, we present a convenient way for rigorously conserving energy. 

2.2 Energy based formulation 

While most computational techniques developed so far only conserve energy in an ap- 
proximate manner [U [9], here we show that an energy-based formulation provides a 
convenient and rigorous way to conserve energy in the relaxation time approximation. 

Adopting a similar approach as in [2j to derive the Equation of Phonon Radiative 
Transfer, we multiply by hw to obtain 

de Joc - 

dt 



+ V g Ve = (10) 
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which we will refer to as the energy-based BTE. Here, e = huf and e 



loc 



loc 



Equation ( 10 ) can be simulated by writing 



8tt s £, 



eff 



E* 3 



(11) 



where £ e tt is defined as the effective energy carried by each computational particle. 



Statement (11) defines computational particles that all represent the same amount of 
energy. From the point of view of phonons, comparing ^ and (11) shows that the 



effective number of phonons represented by the newly defined particles is variable and 
is linked to the effective energy by the relation £ e ff = N e ffhuj. By analogy with 



the description of section 2.1, computational particles defined by (11) obey the same 



computational rules as in the previous Monte Carlo approaches. Modifications appear 
at three levels: 

• When drawing particle frequencies during initialization, emission from boundaries 
or scattering, the distribution functions that we use must account for the factor 
hu. For example, when initializing an equilibrium population of particles at a 
temperature T , one has to draw the frequencies from the distribution 



cxp 



k b T 



1 



(12) 



• Calculating the energy in a cell is straightforward and simply consists in counting 
the number of computational particles. The energy associated with N particles is 
given by £ e ffN. 

• Since the energy in a cell is proportional to the number of particles, there is no need 
for an addition/deletion process: energy is strictly and automatically conserved by 
simply conserving the number of particles. 



2.3 Deviational formulation 

In this section we introduce an additional modification which dramatically decreases the 



statistical uncertainty associated with Monte Carlo simulations of (10). Our approach 
belongs to a more general class of control- variate variance reduction methods for solving 
kinetic equations [FT] [T9] in which the moments < R > of a given distribution / are 
computed by writing 

J Rfdxdc = J R(f - f eq )dxdc + J Rf eq dxdc (13) 

where the first term of the right hand side is computed stochastically and the second 
term is computed deterministically. If f eq ~ /, the variance reduction is large because 
only a small term is determined stochastically (see Figures [l]and[2]). 
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In the present context, this methodology provides significant computational savings 
when an equilibrium (constant temperature) state exists nearby, which is precisely the 
regime in which statistical noise becomes problematic (low signals). The degree of 
variance reduction achieved by this method is quantified in section [5] 

Let 

4" M = (14) 

" «P(*)"1 

where T eq ^ Ie<j(x, t). Then, it is straightforward to show that e d = e — is governed 
by 

8e d A (e l ° c -ex )~ ^ 

-k +Y * Ve = t (15) 

Therefore, by analogy to the standard particle methods for solving the Boltzmann 
equation, we define computational particles by: 

e d = e - e e T q eq « Svr 3 ^ s(i)5 3 (x - Xi )<5 3 (k - ki)6 PJH , s{i) = ±1 (16) 



We will refer to these newly defined computational particles as deviational particles. 



Clearly, deviational particles may be negative since e — e^? can be a negative quantity. 



This is accounted for in the sign term in equation (16). In what follows, we derive 
evolution rules for deviational particles based on (15). 




Figure 1: In standard particle methods, the moments of the distribution are stochasti- 
cally integrated. 



3 Algorithm 

The variance-reduced algorithm is very similar to its non-variance reduced counterpart 
and comprises an initialization step followed by a splitting algorithm for time integra- 
tion. The main change lays in the distributions from which deviational particles are 
sampled. 
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Figure 2: In a control- variate formulation, the stochastic part is reduced to the calcu- 
lation of the deviation from a known state, which is much smaller. 



3.1 Initialization 

The algorithm proceeds by choosing the equilibrium state at temperature T eq from which 
deviations will be simulated. Although this choice can be quite critical in the efficiency 
of the method (the smaller the deviation from the chosen equilibrium state, the smaller 
the number of deviational particles required for a given statistical uncertainty, or for a 
fixed number of deviational particles, the larger the variance reduction), it is usually a 
natural and intuitive choice. 

In some cases, the equilibrium state is the same as the initial state. In such a situ- 
ation, the simulation starts with no particles. Nevertheless, one still has to choose the 
deviational effective energy for subsequent use. In the various examples discussed 
below, this parameter was chosen as follows: based on a guess of the upper bound on 
the deviation of temperature at steady state, the deviational energy of the system can 
be estimated using 



AE 



ui=0 



1 



1 



exp 



ftw 
k b T 



1 



dui 



(17) 



This estimate of the deviational energy allows £^ff to be (approximately) determined 
based on the desired number of computational particles. 

If the initial state / is different from the equilibrium distribution, particles need 
to be initialized in the computational domain. Their frequencies and polarizations are 
drawn from the distribution 



D(uj,p)e d (uj) = HujD(u),p) 



f° 



1 



(18) 



Typically, f° is an equilibrium distribution at some temperature T, whereby the above 
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expression reduces to 



D{uj,p)e d (u) = huD(uj,p) 



exp {fir) 



(19) 



This function is positive if T > T e<? and negative if T < T e(J . As a result, in the latter 
case, particles are assigned a negative sign. Drawing the frequencies is performed as 
in [8], namely by subdividing the frequency range in bins (generally, about 1000 bins 
are considered enough), defining a discretized and normalized cumulative distribution 



from (19), uniformly drawing a random number between and 1 and finding the bins 



it corresponds to in order to match the normalized cumulative distribution. 



3.2 Advection 



Since the left hand side of (15) is analogous to that of ([!]), the advection substep is 
unchanged. In other words, during the time step At, particles of group velocity V g {uj,p) 
are simply advected by ~V g (uj,p)At. 



3.3 Sampling substep 

Sampling the local temperature and pseudo-temperature requires a few changes from 
the non-variance reduced method, namely 

• Let Cj be the set of indexes corresponding to the particles inside cell j of volume Vj 
at time t. Since each particle represents the same amount of energy, the deviational 
energy is given by 



(20) 



where J\fj~ and J\fj are respectively the number of positive and negative particles 
inside the cell j. 

• The corresponding temperature Tj is then calculated by numerically inverting the 
expression 



£J=0 



1 



exp 



1 exp 



( ^ 

V k b T el 



duj (21) 



Similarly, once Tj is known, the deviational pseudo-energy is computed using 



T(Ui,Pi,Tj) 



(22) 
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The corresponding pseudo-temperature [T/ oc ] j is calculated by numerically invert- 
ing 



V-i 



D(uj,p)huj 



1 



do; (23) 



3.4 Scattering step 

During the scattering step we integrate 



dt 



T(u,p,Tj) 



for a timestep At, where 



„loc 



l l 

.^(to)- 1 " 6 ^^)- 1 . 



(24) 



(25) 



We select the particles to be scattered according to the scattering probability (specific 
to each particle's frequency and polarization, and depending on the local temperature) 



P(ui,Pi,Tj) = 1 - exp ( At ) 



(26) 



The pool of selected particles represent a certain amount of deviational energy £ ^ (JV^j — 
7V S ~ ), where A/" s ~j and A/"~j refer respectively to the number of positive and negative se- 
lected (i.e. scattered) particles in cell j. This pool of selected particles must be replaced 
by particles with properties drawn from the distribution 



D(uj,p)(e 



loc 



-, c <j 



T(u!,p,Tj) 



D{u,p)huj 



1 



1 



1 



(27) 



which is either positive for all frequencies and polarizations or negative for all frequencies 
and polarizations. In other words, scattered particles must be replaced by particles 



which all have the same sign as e — and which respect the energy conservation 
requirement. Therefore, out of the A/" s + j +A/"~,- selected particles, we redraw properties 



for 



8,J 



AC 



particles. The 



and are given the sign of e 



of them according to the distribution ( 27 ) and delete the other selected 



A/?, -AT 



particles to be kept are chosen randomly inside the cell j 



J.OC 



This process tends to reduce the number of particles in the system and counteracts 
sources of particle creation within the algorithm (e.g. see boundary conditions discussed 
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in the next section) . A bounded number of particles is essential to the method stability 
and the reduction process just described is a major contributor to the latter [TTJ [12]. 
Hence, in a typical problem starting from an equilibrium state that is also chosen as the 
control, the number of particles will first increase from zero and, at steady state, reach 
a constant value that can be estimated by appropriately choosing 6^ as described in 
section 



3.1 The constant value will usually be higher (but of the same order) than the 



estimated value: indeed, the rate of elimination of pairs of particles of opposite signs 
depends on the number of particles per cell and therefore on the spatial discretization 
chosen (the finer the discretization, the smaller the number of particles per cell and 
therefore the smaller the rate of elimination). 



3.5 Boundary conditions 

In phonon transport problems, various types of boundary conditions appear. Isothermal 
boundary conditions, similar by nature to a black body, have been used in several studies 
[5J, [9] . Adiabatic boundaries also naturally appear [HI [20]. Recently, a class of periodic 
boundary conditions has also been introduced [10]. The deviational formulation adapts 
remarkably well to these different classes of boundary conditions. 



3.5.1 Adiabatic boundaries 

Adiabatic boundaries reflect all incident phonons. This reflection process can be divided 
into two main categories: diffuse reflection and specular reflection. In both cases, it 
is assumed that the polarization and frequency remains the same when a phonon is 
reflected. The only modified parameter during the process is the traveling direction. 

i Specular reflection on a boundary dn of normal vector n can be expressed, in 
terms of energy distribution, by 

e(x,k) = e(x,k') (28) 

where k' = k — 2(k • n)n and x 6 dn. Since the equilibrium distribution e^? is 
isotropic, then substracting it from both sides simply yields 

e d (x,k) = e d (x,k') (29) 

In other words, deviational particles are specularly reflected 

ii Diffuse reflection amounts to randomizing the traveling direction of a phonon 
incident on the boundary, in order for the population of phonons leaving the 
boundary to be isotropic. Since an equilibrium distribution is already isotropic, 
incident deviational particles are treated identically to real phonons. 
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3.5.2 Isothermal boundaries 



In the case of an isothermal boundary at temperature T&, incident phonons are absorbed, 
while the boundary itself, at temperature T&, emits new phonons from the equilibrium 
distribution corresponding to T5. The emitted heat flux per unit radial frequency is 
expressed by 

// 1 D{u,p)V g (uJ,p)hui 

~ 4 ^ expf^-1 { ) 

Substracting the heat flux per unit radial frequency corresponding to a boundary at 
equilibrium temperature, we obtain 

Ct = \ E D ( ., V) V M ( - 7 ^ T - ) (31) 



which gives the frequency distribution of emitted particles. Traveling directions must 
be chosen accordingly, as explained for example in [8]. 

3.5.3 Periodic unit cell boundary conditions 

Heat transfer in periodic nanostructures is a subject of considerable interest in the 
context of many applications. Such nanostructures are considered in Hao et al. [TP] , 
in Huang et al. |21j and in Jeng et al. [4j. Hao et al. developed periodic boundary 
conditions that allow efficient simulation of such structures by considering only one unit 
cell (period). In this section we review the work of Hao et al. [10J and explain how 
the deviational particle formulation presented here lends itself naturally to this type 
of boundary condition. Simulations using these boundary conditions are presented in 
section 16.11 

We consider a 2D periodic structure depicted in Figure [3] in which square unit cells 
containing two rectangular voids are organised in a square lattice. Our interest focuses 
on determining the effective thermal conductivity of such a structure as a function of 
d, the degree of alignment. 

The formulation introduced by Hao et al. amounts to stating that, at the boundaries, 
the deviation of the phonon distribution from the local equilibrium is periodic. Using 
the notations from Figure [3j this condition can be written as 



f + _ f e i — f+ _ f e i 
f — — f e i — f - — f e< 

J l JTi ~ J 2 JT 2 



T ? (32) 



where f^? and f^? refer to the equilibrium distributions at temperatures T\ and T2, 
where the superscript + denotes particles moving to the right (with respect to figure [3]) 
and where superscript — refers to particles moving to the left. This formulation enforces 
at the same time the periodicity of the heat flux and a temperature gradient. In terms 



12 





■l 


■l 


'l 






■l 


■l 


■l 


L 




■l 


■l 


'l 



Periodic 
flux 



fl 



/ \ 


1 






1 



Jt 



Figure 3: Example of a periodic nanostructure. Each periodic cell comprises two rect- 
angular voids with diffusely reflecting walls. This nanostructure, and in particular the 



influence of the parameter d, is studied further in section 6.1 



of deviational energy distributions, this relation becomes 

f eq 



M/i + - fZ 



rz = m/ 2 + - fz - %) 



M/r - rz - fZ = m/ 2 - - fZ - fZ 



(33) 



which amounts to 




e 2 



d,+ eg 
e T 2 



d,-_ eg (34) 

Computationally, this formulation can be implemented by emitting new particles 
from both sides while periodically advecting the existing particles. Without any loss of 
generality, let us assume that T\ > T 2 . Particles emitted from the hot side originate 
from the distribution 



Therefore, at a given point on the boundary, denoting 9 the angle with respect to the 
normal and 4> the azimuthal angle, the flux per unit radial frequency locally emitted 
from boundary 1 ("hot" side) in the solid angle d£l = sin 0d9d(j) can be expressed as 



it 



p 

.D(u,p) 



47T 



-V g {ijj,p) cos 9 sin 



E 



4tt 



V g (uj,p) cos 9 sin 



+ (4! 



D eg 
-T 2 > 



DM) 

4-7T 



Vg(u,p) cos 9 sin 9d9d(f> 



crossing boundary 2 



new particles generated 



(36) 

Similarly, the flux per unit radial frequency locally emitted from boundary 2 ("cold" 
boundary) can be expressed as 



v v — 



d-D(uj,p) . 
e/ — : Vg(uj,p) cos 9 sm 



4tt 



(e eg 



DM) 



4tt 



V g (oj, p) cos 9 sin 9d9d(j) 



crossing boundary 1 



new particles generated 



(37) 
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Hence the boundary condition can be enforced by: 

i Moving all particles and applying periodic boundary conditions to those crossing 
a periodic boundary: a particle leaving the system on one side is reinserted on the 
other side. 

ii Generating new particles from the distribution 

(4!-4*)^^p) (38) 



The number of new particles is given by integrating (38) over all frequencies and 
polarizations and by multiplying the result by ir to account for the integration over 
the solid angle f^ Q \q=q cos # sin OdOdcfi. The traveling direction of these particles 
is randomized on the half-sphere pointing into the domain and in the case of the 
hot boundary they are sent traveling to the right with a positive sign. Taking 
their mirror image, negative particles with the same properties are emitted by the 
cold boundary. 

4 Validation 

4.1 A ballistic problem 

In order to validate the proposed formulation, we first consider a one-dimensional system 



bounded by two isothermal (3.5.2) boundaries that are sufficiently close - their distance, 
L, is much smaller than all phonon mean free paths - that transport can be modeled 
as ballistic. The system is initially at a uniform equilibrium temperature To, when at 
t = + the temperature of the isothermal walls impulsively changes to To =L AT. 

Appendix [B] presents an analytical solution for the resulting transient evolution of 
the temperature field that is used here for comparison with our simulations. A particu- 
larly interesting case is the Debye model which, when coupled with small temperature 



amplitudes, allows a linearization of the general relation ( 55 ) to provide a fairly simple 
closed- form solution (56). Figure [4] shows a comparison between this solution and the 
variance-reduced Monte Carlo result. The simulation was run with T eq = Tq and the 
phonon velocity was taken to be 12, 360m. s -1 [10]. Excellent agreement is observed. 

4.2 Heat flux and thermal conductivity in a thin slab 

In this section we continue to validate our formulation by calculating the thermal con- 
ductivity of a thin silicon slab bounded by two diffusely reflecting walls a distance d 
apart in the z direction (see Figure [5]). The slab is infinite in the x and y directions. 

This problem is considered here because the solution can be expressed analytically. 
We introduce the local deviation function f d = f — f loc and, denoting the temperature 
gradient by dT/dy, rewrite the BTE at steady state as: 

Afloc jrp rd 
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Figure 4: Transient temperature profile in a one-dimensional ballistic system whose 
boundary temperatures undergo an impulsive change at t = 0. Initially, the system is 
in equilibrium at temperature To = 300K. At t = + , the wall temperatures become 
T ± AT; here, AT = 3 K. 




Figure 5: Heat conduction in a silicon slab due to an imposed temperature gradient in 
the y direction. Slab is infinite in the x and y directions. 
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Figure 6: Simulation geometry. Boundaries at z = and z = L are diffusely reflecting. 
Infinite domain in y direction is terminated using periodic boundaries L = lOOnm apart. 



This equation can be solved to yield, in the coordinate system introduced in Figure 



f d (z,uj,p,9,0 < (jxir) 



A(w,p,T )cos(6>) 



df l ° c {u,T Q ) dT 
dT dy 



f d (z,u,p,9,-ir < cp< 0) 



A(u,p,T )cos(6) 



df loc (u,T ) dT 
dT dy 



1 — exp 



1 — exp 



A(cj, p, To) sin(#) sin(0) 
z — d 

A(cj, p, 7q) sin(#) sin(0) 



(40) 



(41) 



where A(uj,p, To) is the average mean free path at frequency lo, polarization p and 
temperature To, given by 



A(u],p,T ) = V 9 (u;,p)t(lj,p,T ) 



(42) 



Moments of this solution can be numerically integrated to yield values for the heat flux 
and the thermal conductivity of the slab. 

In the simulation, we calculate the thermal conductivity by measuring the steady 
state heat flux in response to a temperature gradient along the y axis (see Figure [6]) . 
Due to the translational symmetry of the system, we impose the temperature gradient 



using the periodic unit-cell formulation presented in section 3.5.3, which allows us to 
use a finite system size in the y-direction, taken to be L = lOOnm. In order to measure 
the thermal conductivity at To, a temperature gradient is imposed by setting a target 
temperature of To + AT for the hotter of the two boundaries and To — AT for the colder 
boundary, and we proceed as explained in 3.5.3 The deviational method allows the 



solution of this problem for AT <C To (here, AT = 0.05K), in contrast to non- variance- 
reduced methods that would require AT ~ Tq to achieve statistically significant results. 
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The best choice for the equilibrium (control) temperature is clearly T eq = Tq = 300K. 
Initializing the simulation at equilibrium at Tq is also convenient, because no particles 
need to be generated for the initial configuration. 

Figure [7] compares the heat flux in the y direction inside a slab of silicon (see Ap- 
pendix A for material parameters) of thickness d =100nm, as computed by the devia- 
tional method, to the analytical solution. Figure [8] compares the thermal conductivity 
of the slab at To=300K as a function of d computed from the deviational method and 
from the analytical expression. Very good agreement is observed in all cases. 




Figure 7: Spatial variation of the axial (in the y direction) heat flux in a thin film 
with a thickness d = lOOnm, computed theoretically and compared to the result of the 
deviational simulation. 



5 Computational efficiency 

The variance-reduced method developed here allows substantial improvement in the 
relative statistical uncertainty, a /AT, compared to non- variance-reduced simulations. 
Here, a is the standard deviation in the temperature measurement and AT is the 
characteristic temperature difference (as, for example, in the validation case studied in 

Figure [9] compares the relative statistical uncertainty of the variance-reduced with 
the standard method. The reported data was obtained by simulating equilibrium at 
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Figure 8: Theoretical values of the thin film thermal conductivity at Tq = 300K, com- 



puted by numerical integration of the theoretical expressions (40) and (41). Comparison 
with the values obtained from the deviational simulation. 
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some temperature T±, and defining AT = T\ — Tq as the characteristic signal that 
needs to be resolved. By choosing T eq = Tq in the deviational method, we ensure finite 
deviation from equilibrium is considered and thus the statistical uncertainty is non-zero. 
Simulating an equilibrium state is a matter of convenience, because in non-equilibrium 
problems the number of particles and thus the local statistical uncertainty varies as a 
function of space in the deviational simulation and is thus difficult to quantify precisely; 
simulations of simple problems (e.g. Couette-type problems) in the past E21 E3] have 
yielded very similar results. We also note that even though a/ AT is strictly speaking 
the ratio of statistical uncertainties, it serves as a good approximation to the ratio of 
computational cost, because the cost of the deviational simulation is very similar to that 
of standard Monte Carlo methods. Specifically, the speedup provided by the deviational 
method is given by the square of the relative statistical uncertainties. 

A very interesting feature of variance-reduced methods is that the standard devia- 
tion of the results is proportional to the amplitude AT of the signal, as shown in Figure 
[9] (see also [22l [23j). As a consequence, variance-reduced methods are able to pro- 
vide the desired relative statistical uncertainty (signal to noise ratio) for arbitrarily low 
signals without requiring more computational effort. In contrast, in the case of the non- 
variance-reduced method, it is more computationally expensive to obtain the desired 
level of relative statistical uncertainty for small variations in temperature, than for large 
variations in temperature. In these methods, for AT « Tq, the statistical uncertainty 
is approximately constant (set by equilibrium fluctuations) and thus a /AT ~ 1/AT. 
As a result, the speedup offered by the variance-reduced methods scales as 1/(AT) 2 . 
For example at AT/T « 1(T 2 (i.e. AT « 3K at room temperature) the speedup is 
approximately 4 orders of magnitude (see Figure [9]); at AT/Tq ss 10 -3 , the speedup is 
approximately 6 orders of magnitude. 

6 Applications 

In this section we present some applications of the deviational method to problems of 
current engineering interest. Modeling work in these areas is still ongoing; the objective 
of this discussion is mainly to showcase the capabilities of the proposed method. 

6.1 Thermal conductivity of nanoporous silicon: influence of nanopore 
alignment 

Decreasing thermal conductivity as a means of improving the thermoelectric effect has 
received considerable attention, and nanostructures are a novel approach towards this 
goal. Similarly to Huang et al. [21] and Jeng et al. [I], we assess here the thermal 
conductivity of novel nanostructured materials. The nanostructure considered here is 
made of rectangular pores as shown in Figure [3j We model it as a 2D problem (possible 
if the material boundaries in the directions normal to the plane shown in the figure can 
be approximated as specularly reflecting). Figure [To| shows the periodic cell considered 
and defines the parameter d that we use to describe the spatial distribution of the 
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Figure 9: Comparison of relative statistical uncertainties for equilibrium systems at 
temperature T\ with AT = T\ — Tq and Tq = 300 A". 



pores. The thermal conductivity in the y direction is measured by imposing periodic 



unit-cell boundary conditions as explained in section 3.5.3 with a temperature difference 
of 0.1K across the unit cell. Using the data of Appendix A, the contributions of the 
different mean free paths to the bulk thermal conductivity can be calculated. A plot of 
the effective thermal conductivity, as computed with the deviational variance-reduced 



method, is displayed in Figure 10 The thermal conductivity is reduced by almost a 
factor of 2 because of this geometrical effect. This highlights the importance of ballistic 
effects. 



The importance of ballistic effects is further highlighted by Figure 11 which shows 
that at To = 300K, mean free paths from 50nm to 10//m contribute significantly to 
the thermal conductivity of the bulk material; the presence of voids with period of 
lOOnm affects the contribution of all mean free paths, but completely suppresses the 
contribution of all mean free paths greater than about one micrometer. Tuning the 
alignment parameter, decreases further the contribution of the mean free paths between 
50nm and 1/xm. 

6.2 Simulation of thermal conductivity spectroscopy 



Figure 12 depicts an experimental setup developed in the MIT Nanoengineering Lab 
[25] as a prototype "thermal conductivity spectroscopy" system. This experiment is 
based on pump- probe transient thermoreflectance, in which a pump pulse is used to 
change the physical properties of a sample and a probe pulse is used to measure the 
change. In this experiment, a thin film of aluminum (thickness between 50 and lOOnm) 
is deposited on a silicon wafer and is initially at uniform temperature, say 300K. At 
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Figure 10: (a) Temperature field in a unit cell of a periodic nanoporous material, (b) 
Thermal conductivity clS ct function of parameter d. 
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Figure 11: (a) Thermal conductivity accumulation as a function of the mean free path. 
The bulk conductivity was computed by numerical integration of the thermal conduc- 
tivity per unit frequency Tv 2 C UJ j?> [EJIMIII]; here, is the heat capacity per frequency 
unit, (b) Normalized thermal conductivity accumulation, highlighting the influence of 
ballistic effects on the thermal conductivity. 
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t = 0, localized laser irradiation creates a hot spot, shown in figure 12 as centered on the 
origin (r = 0, z = 0) of the coordinate system. A reliable description of the subsequent 
evolution of the temperature field is central to interpreting the experimental results and 
creating a means for inferring phonon mean free paths (the goal of this experiment) 
from experimental measurements (e.g. surface temperature) 

Given the scale of the aluminum slab, the impulsive nature of the heating, and the 
short duration of the phenomenon, phonon ballistic behavior needs to be accounted for, 
necessitating a Boltzmann treatment. However, this problem is very difficult (if not 
impossible) to simulate using standard Monte Carlo methods: the initial perturbation 
to the temperature field is small in amplitude (see below) which makes resolution of 
transient results very costly. Moreover, the need to simulate early as well as late times 
and avoid artifacts from artificial domain termination makes the simulation of a large 
computational domain necessary, even though the original hot spot is very small. In 
traditional Monte Carlo methods, this large computational domain would need to be 
filled with particles. 

The method proposed makes this calculation possible. Simulating the deviation 
from equilibrium allows the calculation to proceed using zero particles in regions not yet 
affected by the heating pulse. Thus, in addition to variance reduction which removes 
the limitations associated with statistical uncertainty, simulating the deviation from 
equilibrium simultaneously considerably reduces the computational cost resulting from 
the multiscale nature of this problem. We also note that by taking the equilibrium 
distribution at 300K, the simulation only has positive particles. Hence there will be 
no cancellation of particles and the entire simulation will run with a fixed amount of 
particles. 

In practice, one can exploit the cylindrical symmetry in order to reduce the problem 
dimensionality: the resulting temperature field is expected to depend only on the depth 
z and on the distance from the center of the pulse, r. Therefore, we can use toroidal 
cells to sample the temperature and process the scattering. The only drawback is that 
cells near the center, at small radius, will have a smaller volume and will sample the 
temperature over a smaller number of particles, thus yielding noisier results in these 
regions. 

6.2.1 Initial condition 

As stated above, since the material is originally at equilibrium at To = 300K it is most 
convenient, but also computationally efficient, to choose T eq = Tq. Laser irradiation 
introduces a heating effect in a thin layer close to the irradiated surface which has been 
parametrized [15J using the following expression 

2r 2 



AT(r,z) = T-T = Ae W [--2-f3z) (43) 



with A = IK, i?o = 15/um and j3 1 = 7nm. This expression is used here as an initial 
condition for the material temperature. Regions for which AT < .005-fT were taken to 
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be at equilibrium at To = T eq (no particles). 
6.2.2 Interface modeling 

The top surface of the aluminum material (z = 0) is modeled as a diffusely reflecting 
wall. 

Modeling the interface between the two materials accurately is still an active area 
of research. Here, we chose to use a recently developed model [15, 26J which relates the 
transmissivity to the interface conductance G through the expression 

< Pi^CtV^ >= i 2 1 (44) 

<C 1 y 9 ,i> + <C 2 V g , 2 > + 2G 

Here, Pi-^j denotes the probability for a phonon to pass through the interface from 
material i to j; the brackets denote integration over frequency and sum over polarization, 
while C\ and C2 denote the volume heat capacity per unit frequency in media 1 and 2, 
respectively. In this model, we assume that the interface is totally diffuse: the direction 
of an incident particle is reset regardless of the transmission or reflection of the particle, 
while its frequency and polarization are retained |27j . For the interface conductance G, 
we use the experimental value G = 1.1 x 10 8 Wm _2 K _1 [26 . 
We also utilize the expression [27] 

D 1 (u,p)V Xt g(u, P )f^P 1 ^2{u},p) = Ih(u,p)V 2 ,g(u,p)f%P^x(u,p) (45) 

which relates the probability for a phonon with radial frequency u and polarization p 
to pass through the interface from 1 to 2 to the probability to pass from 2 to 1. 



We can easily verify that relation ( 45 ) applies when the deviational energy e is used 



instead of the phonon distribution. Additionally, expression (44) which relies, among 



other things, on (45) [151 [26], also remains unchanged when applied to deviational 
particles. 

Following \15\ I26j. we let Pi_>2 be a constant (which makes it easy to calculate from 



(44)) and deduce J*j-n from (45). In our case we chose to set Pai-^Si constant, except 
for the high frequency transverse acoustic modes; since the cutoff frequency of the TA 
branch in Si is lower than the TA cutoff frequency in Al, phonons with such frequencies 
must undergo total reflection [15J. Similarly, LA phonons in Si whose frequency is above 
the aluminum LA branch cutoff frequency are totally reflected. 



6.2.3 Domain termination 

At long times, phonons may travel far from the hot spot. In order to avoid discretizing an 
infinite domain with computational cells (for calculating the temperature) we restrict our 
discretization to a finite (but large) "nominal" domain. In order to simulate accurately 
and consistently the actual system, we keep track of the particles even after they have 
left the nominal part of the domain. 
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Particles that leave this domain are not sampled (for calculating the temperature 
and pseudo temperature), but are still scattered by assuming a local temperature of 
300K as an input parameter for the relaxation time. This amounts to a linearization 
of the collision operator at T = 300K and is based on the reasonable assumption that 
sufficiently far from the heating source, the temperature is very close to 300K. Particles 
that leave the nominal part of the domain may reenter the nominal domain, hence 
ensuring a rigorous treatment of the semi-infinite region. 

Particular care is taken to ensure that the frequency and polarization of a particle 
is drawn from the correct distribution, because energy conservation — built into the 
simulation method — requires that the number of particles is conserved by the scattering 
process and is inconsistent with approximations which do not conserve energy. For 
example, setting T\ oc =300K is inconsistent with energy conservation because e loc (Ti oc = 
300K) — e eq = 0, which implies no particle generation, which in the presence of particle 
deletion due to the term —e d /r leads to net particle and thus energy loss. This situation 
can be rectified by allowing the temperature at the particle position to be different from 
T eq ; specifically, we write T = T eq + e and expand 

T(u,p,T eq ) r(uj,p,T eq ) dT 

Frequencies and polarizations are thus drawn from 



D(co,p) 04? 



r(uj,p,T eq ) dT 



(47) 



since (46) once normalized, does not depend on the local e. As before, energy conser- 
vation is ensured by simply conserving the particles. 

In addition to providing a method for terminating simulations, this approach rep- 
resents a promising avenue for treating the entire simulation domain in the limit that 
linearization of the collision operator is appropriate. The advantage of this formulation 
is significant reduction in computational cost because evaluation of the local temper- 
ature and pseudo-temperature is not required every timestep. Further details will be 
presented in a future publication. 

6.2.4 Simulation results 



Figure 13 and 14 show that the variance-reduced method developed here can calcu- 
late the temperature field with small statistical uncertainty. This is remarkable given 
the minute temperature differences (O(0.01)K) present in this problem, especially at 
late times. For such temperatures, according to Figure [9| the speedup compared to a 
standard Monte Carlo method is on the order of 10 9 . 



Figure 14 compares our simulation results with a numerical solution of the heat 
conduction equation (Fourier's Law). The differences between the two predictions are a 
result of non-diffusive (ballistic/transitional) effects. The detailed information available 
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in simulations of this phenomenon can assist in the development of methodologies for 
characterizing carrier mean free paths from comparisons such as the one shown in Figure 
14. Here, we note that the present calculation does not account for thermal transport 
by electrons in aluminum. This was neglected in the interest of simplicity and because 
the primary focus of this experiment is transport through the silicon substrate |15j . 
Thermal transport by electrons in aluminum will be considered and evaluated in a 
future publication. 




Figure 12: System composed of a slab of aluminum on a semi-infinite silicon wafer, used 
for transient thermoreflectance (TTR) experiments. At t=0, a laser pulse induces a 
temperature field T(r, z, t). The temperature field evolution after the pulse is computed 
by assuming that the aluminum surface is adiabatic. 
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Figure 13: Variance-reduced temperature field in aluminum slab and the silicon wafer 
after initial heating by a laser pulse. The picture shows the aluminum slab (lOOnm 
thickness) and a portion of the silicon wafer (lOOnm thickness). 
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Figure 14: Surface temperature at the hot spot (averaged over the region < r < 2/xm, 
< z < 5nm) as a function of time after initial heating by laser pulse. The difference 
with the solution based on the Fourier model is a result of ballistic effects. 
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7 Discussion 



We have shown that efficient and accurate algorithms for solving the BTE with signif- 
icantly reduced statistical uncertainty can be developed by focusing on the deviation 
from a nearby equilibrium within an energy-based formulation. The energy-based for- 
mulation facilitates exact energy conservation thus improving the simulation fidelity, 
while the variance reduction is made possible by the deterministic information inherent 
in the Bose-Einstein distribution which describes the nearby equilibrium. The proposed 
method was validated using analytical solutions of the Boltzmann Transport Equation. 
Very good agreement with the analytical results was found. 

The proposed algorithm was used to study the effect of porosity on the effective 
thermal conductivity of pure silicon. Our results show that staggering periodically 
arranged voids at small scales exploits ballistic shading to effect reduction in the effective 
thermal conductivity. A more systematic investigation of the effects of porosity on the 
effective conductivity of silicon — including anisotropic effects — will be the subject of 
future work. 

We also presented simulations of a recently developed experimental technique known 
as thermal conductivity spectroscopy, in which the transient response of a thin aluminum 
slab over a silicon wafer to a localized heating induced by a laser pulse is used to infer 
properties of heat carriers. This simulation required the development of a domain 
termination algorithm for rigorously treating deviational particles as they travel to 
regions far from the heating source, without having to sample these particles everywhere 
in this semi-infinite region. The algorithm developed corresponds to a linearization of 
the collision operator and may, in fact, form the basis of a significantly more efficient 
simulation approach valid in cases where linearization is appropriate. 

In addition to illustrating the benefits of variance reduction, simulations of the ther- 
mal conductivity spectroscopy problem also showcase the value of the proposed simu- 
lation approach as a new multiscale method: in contrast to typical multiscale methods 
which focus on spatially decomposing the domain into the particle and continuum sub- 
domains, the present algorithm achieves a seamless transition from one description to 
the other by instead algebraically decomposing the distribution function into a part 
described by particles and a part described deterministically |28j. Although here the 
simplest such implementation has been presented (deterministic description is equi- 
librium at temperature To ^ 7o(x, t)), deviational algorithms featuring a deterministic 
description that varies as a function of space (e eq = e e<? (x)) have been developed [T2"|[T3] 
and shown to achieve improved variance reduction as Kn — > [13j . albeit at the cost of 
a moderately more complex algorithm. In the problem considered here, the continuum 
behavior at large distances from the heat source is in fact equilibrium at To and thus 
the present algorithm is sufficient. However, in other problems where a local equilib- 
rium is present in large parts of the domain, algebraic decomposition using e eq = e eq (x) 
will be able to provide considerable computational savings by considerably reducing the 
number of particles required for its simulation. 
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A Numerical data for scattering rates 



In the simulations presented here we use data for the dispersion relations and for the 
relaxation times of phonons in Al and Si. Dispersion relations are adapted from the 
experimentally measured dispersion in the [100] direction ([29] for Al, [30|. 115] for Si). 

For Al, as in [26, ,15], we assume a constant relaxation time chosen to match the 
desired lattice thermal conductivity. We therefore take the value 

T Al = 10~ n s (48) 

For Si, we use the expressions from [3T], with constants from |15| . Relaxation times 
for acoustic modes are given by 



phonon-phonon scattering 


, LA 


rl 1 


= A l uj 2 T 1A9 exp 


(4) 


phonon-phonon scattering 


, TA 




= A T uj 2 T im exp 




impurity scattering 




rj 1 


= A i oj* 




boundary scattering 




rs 1 


= w b 





where the constants take the following values 



Parameter Al At Aj Wb 

Value (in SI units) 2 x 10" 19 1.2 x 10~ 19 80 3 x 10~ 45 1.2 x 10 6 



The total relaxation time for a given polarization is obtained using the Matthiessen rule 

r-^XVf 1 (49) 

i 



Optical phonons in Si are considered immobile (Einstein model). Einstein's model 
states that the contribution of optical phonons to the vibrational energy per unit volume 
in a crystal is given by pQ 

V[eM^E/k b T) - 1] 1 ' 

where N p = 3 is the number of polarizations, N' = 1 is the number of optical states per 
lattice point, is the Einstein radial frequency (oje = 9.1 x 10 13 s _1 [301 US]), V is the 
volume of a lattice point (with a lattice constant a = 5.43A, V = a 3 /4 = 4 x 10 _29 m 3 ). 
For the relaxation time of optical phonons, we use the value [32J 

t = 3 x 10~ 12 s (51) 
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B Derivation of the transient ballistic ID solution 



Following the impulsive change of temperature at the walls from To to 2] = To + AT 
and T r = Tq — AT, thermalized phonons at temperature T r and 7] are emitted from the 



"right" and "left" wall, respectively (see Figure 15). 



For some arbitrary location x, for a given frequency, polarization and time, the angu- 
lar space can be divided into 3 distinct domains characterized by two angles 9 r (x,ui,p, t) 
and Oi(x, u),p, t) as depicted in Figure 15 Phonons described by < 9 < 9[ were emitted 
by the left wall at a time t > 0. Phonons described by 9\ < 9 < tt — 9 r have been present 
in the system since t=0. Phonons described by tt — 9 r < 9 < tt were emitted by the 
right wall at a time t > 0. 
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Figure 15: At a given point in space, the solid angle can be divided into three distinct 
regions in which the distribution of phonons is known; here, 9\ is given by cos(f^) = 
x/(V g (uj,p)t), while 9 r is given by cos(6* r ) = (L — x)/(V g (uj,p)t) 
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The energy can therefore be written as 



E v (x,t) = -J2 



(x,w,p,t) 



ui Je=o 



e e rp(uj)D(uj,p) sm(9)d9du>... 



f!T—9 r (x,uJ,p,t) 

II e%?(uj)D(uj,p)sm(6)d6duj... 

loj J6=6 l (x,u>,p,t) 



+ 

From geometrical considerations 

L — x 



U! J 9=ix —6 r (x ,u ,p,t) 



e? r (w)-D(w, p) wa.{6)dBdu \ (52) 



cos (6 r (x,uj,p,t)) = min I 1, 



Vg(u,p)t 



1-1 



cos (6i(x,uj,p, t)) = min I 1, 



V g (u,p)t 



1-1 



L — x 

Vg{u,p)t 



Vg(u),p)t 



H I 



L — x 

V g (u,p)t 



(53) 



H I 



Vg(u,p)t 



(54) 



where H is the Heaviside function. Proceeding to the integration in 9, the energy density 
is given by 
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e e ^ o (oo)D(uj,p)duj , > (55) 



The temperature T = T(x, £) is obtained by numerically finding the Bose-Einstein 
distribution corresponding to this energy density. 

Using the Debye model and considering small temperature changes (\T r — Tq\ « To 
and \Ti —Tq\ « To), the resulting temperature field can be expressed in a simpler form. 
The first assumption allows the removal of the frequency and polarization dependence 
on the group velocity, while the second assumption allows the linearization of the Bose- 
Einstein terms in the integrals. Several simplifications can then be carried out to yield 
the following expression for the temperature field 
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